In [1]:
# Lee Mouse WT scRNA-seq data figures 
# 05-20-2025
# Minhyung Kim
library(Seurat)
library(ggplot2)
library(dplyr)
set.seed(2025)
path <- "~/projects/2025/LeeNaeun/"
path_out <- "~/projects/2025/LeeNaeun/out/"
source("~/projects/MK_Rcode/volcano_scdata.R")

Loading required package: SeuratObject

Loading required package: sp


Attaching package: ‘SeuratObject’


The following objects are masked from ‘package:base’:

    intersect, t



Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
In [2]:
if (!file.exists(paste0(path,"Lee_WT_052025.rds"))) {
  path_1 <- "~/projects/2025/LeeNaeun/HN00233422_result_10X/"
  path_2 <- "~/projects/2025/LeeNaeun/HN00233483_result_10X/"
  B1 <- c(1, 2)
  B2 <- 3
  f_name <- "/count/sample_filtered_feature_bc_matrix"
  Lee_Mouse <- list()
  CITE <- list()
  for (i in B1) {
    f_path <- paste0(path_1, i, "/per_sample_outs/", i, f_name)
    data <- Read10X(data.dir = f_path)
    Lee_Mouse[[paste0("sample_",i)]] <- CreateSeuratObject(counts = data$`Gene Expression`, project = paste0("sample_",i))
    Lee_Mouse[[paste0("sample_",i)]]$batch <- "B1"
    CITE[[paste0("sample_",i)]] <- data$`Antibody Capture`
  }
  for (i in B2) {
    f_path <- paste0(path_2, i, "/per_sample_outs/", i, f_name)
    data <- Read10X(data.dir = f_path)
    Lee_Mouse[[paste0("sample_",i)]] <- CreateSeuratObject(counts = data$`Gene Expression`, project = paste0("sample_",i))
    Lee_Mouse[[paste0("sample_",i)]]$batch <- "B2"
    CITE[[paste0("sample_",i)]] <- data$`Antibody Capture`
  }
  Sobj <- merge(Lee_Mouse[["sample_1"]], y = list(Lee_Mouse[["sample_2"]], Lee_Mouse[["sample_3"]]), add.cell.ids = paste0("sample_",c(1:3)))
  saveRDS(Sobj, file = paste0(path,"Lee_WT_052025.rds"))
  saveRDS(CITE, file = paste0(path_out,"Lee_WT_CITE_052025.rds"))
  rm(path_1, path_2, B1, B2, f_name, Lee_Mouse)
}
if (!exists("Sobj")) {
  Sobj <- readRDS(paste0(path,"Lee_WT_052025.rds"))
  CITE <- readRDS(paste0(path_out,"Lee_WT_CITE_052025.rds"))
}
In [3]:
CITE

  [[ suppressing 34 column names ‘AAACCAAAGACAATTC-1’, ‘AAACCAAAGACGACGA-1’, ‘AAACCAAAGATGTATC-1’ ... ]]

  [[ suppressing 34 column names ‘AAACCAAAGACTATGT-1’, ‘AAACCAAAGCACAATT-1’, ‘AAACCAAAGCGAAGAT-1’ ... ]]

  [[ suppressing 34 column names ‘AAACCAAAGACTAATC-1’, ‘AAACCAAAGATTCGAT-1’, ‘AAACCAAAGCATCTTC-1’ ... ]]


$sample_1
7 x 20092 sparse Matrix of class "dgCMatrix"
                                                                              
CD138     1   .   2   3   1  1   3   4  1   1  1   1  4   2   2  10   4   2  7
CD19    109 177 124 114 187 69 129 162 21  22 13 160 83 235 135 148 324 151 77
CD11b     3   2   4   1   6  4   5  11  .  44  4   2  3   5   3   3   4   6  1
CD11c     8   4   5   1   4  6   3   5  4  16  1   8  4   2   2   3   9  23  3
CD43     14  17  23  14  12 84  38  50 22 105  4   8 25  56  54  59 257  18 28
CD5      26  16  20  18  26 49  31  38 18  18 16  12 18  30  84  20 263  24 12
CD21_35  39 133  88  42  58 28  52  89  8   1  1 191 42  63 137  68  12  55 21
                                                                        
CD138     1   1   6  2   5  13  2   9   3   2    4  1   1   6   5 ......
CD19    124 112 174 89 272 258 93 568 101 199  407 59 174 158 143 ......
CD11b     5   3  18  1   4   5  3   6   3   6    7  2   .   8   2 ......
CD11c     2   1   5  3   9   4  6  35   3   5   38  1  11   3   7 ......
CD43     33  79  23 26  63  59  7  65  16  29   58 10  13  53  23 ......
CD5      20  18  37 17  32  27  9  56  84  18   44 17  32  45  41 ......
CD21_35  63  77  43 31 273 190 34 808  48 179 1104  8 231  87  49 ......

 .....suppressing 20058 columns in show(); maybe adjust options(max.print=, width=)
 ..............................

$sample_2
7 x 19107 sparse Matrix of class "dgCMatrix"
                                                                             
CD138      7   5   2   2   8   5 21   2  1   3   8   4   5   7   4   4   6  4
CD19      22 221 135 146 226 138 80 105 84 137 266 104 100 248 222 190 252 37
CD11b     58 558   3   5   .   4  .   3  1   1   3   4   2   8   3   3   3  1
CD11c   1194  12   1   5   2   2  3   3  1   4   1   .   1  17  12   3   4  2
CD43      75 595  36   8   6  33 33  33  5  39 154  24  14 111  43 206  20  .
CD5      244  73  28  17  23  37 23  54 15  23  49  24  25 202  47  87  38 18
CD21_35   11  77  54 124 181  72 49  75 46  35  94   8  77  57 124   4 306 67
                                                                            
CD138     1  7  12   3   7  3   5   9   2   1   1   4   1  11   5   4 ......
CD19    186 49 117 123 117 83 103 167 250 174 202 188 132 218 121 127 ......
CD11b     4  5   1   2   .  3   5   2   5   4   2   2   2   5   2   2 ......
CD11c     3  3   3   3   2  1   4   1   3  10   5   3   1   1   5   2 ......
CD43      9 18  21  41  10  4  34  21  23  27  15  35  13  21  19  23 ......
CD5      35 29  24  19  15 27  43  14  66  30  48  34  22  53  19  17 ......
CD21_35 159 17  34  59 130 44  48  62  86  67 238 179 101 176  90 105 ......

 .....suppressing 19073 columns in show(); maybe adjust options(max.print=, width=)
 ..............................

$sample_3
7 x 24589 sparse Matrix of class "dgCMatrix"
                                                                             
CD138    1   6  21  13   7   6   2   8  13   6   7  97   5  6   4   .   2  15
CD19    84 152 221 569 176  81 124 183 113 161 100 150 134 14 201 131 155 189
CD11b    7   7   5 133   6   2   6  44  13   5   3   2   8  1   6   2   6  10
CD11c    1   2   5   9   3   1   3   5   2   1   5   3   4  2   7   3   2   7
CD43    13  36  55 522  29   5  15  70  10  14  53  31  19  7  13   4   6  32
CD5     18  14  30 245  30 109  19  40  27  18  94  18  15 13  24  11  43  44
CD21_35  2 125  31  39 122  17  75 106  59  73  15  73  69  4  60  73 367  99
                                                                           
CD138     6  4   7  12    7  3   6  11   2   7   7  12   3   3 13  7 ......
CD19    161 99 137  27   33 16 111 344 105 145  15 254 196 111 11 39 ......
CD11b    13  3   4 294 1926  7   4  15   2   4  32   9   3   2  5  9 ......
CD11c     6  2   3 257   67  2   8   3   3   4   7   2   3   .  5  5 ......
CD43      7 25   8 617 3304  5  30 165   6  10 120  75   1  19  9 13 ......
CD5      81 11  16  52   56 11  26  53  26  15  16  26  17  19  6 28 ......
CD21_35 322 86  38  14   12  7  24 112  32  67   4 201 133  45 12 21 ......

 .....suppressing 24555 columns in show(); maybe adjust options(max.print=, width=)
 ..............................
In [4]:
Sobj <- JoinLayers(Sobj)
dim(CITE$sample_1)
dim(CITE$sample_2)
dim(CITE$sample_3)
  1. 7
  2. 20092
  1. 7
  2. 19107
  1. 7
  2. 24589
In [5]:
bb <- cbind(
  CITE$sample_1,
  CITE$sample_2,
  CITE$sample_3
)
colnames(bb) <- colnames(Sobj)
In [6]:
length(colnames(Sobj))
head(colnames(Sobj))
63788
  1. 'sample_1_AAACCAAAGACAATTC-1'
  2. 'sample_1_AAACCAAAGACGACGA-1'
  3. 'sample_1_AAACCAAAGATGTATC-1'
  4. 'sample_1_AAACCAAAGCCAACCG-1'
  5. 'sample_1_AAACCAAAGCCGATTA-1'
  6. 'sample_1_AAACCAAAGCTGCGTG-1'
In [7]:
Sobj

An object of class Seurat 
33696 features across 63788 samples within 1 assay 
Active assay: RNA (33696 features, 0 variable features)
 1 layer present: counts
In [8]:
# QC 
Sobj[["percent.mt"]] <- PercentageFeatureSet(Sobj, pattern = "^mt-")
In [9]:
VlnPlot(Sobj, features = c("nFeature_RNA"), pt.size = 0.1)
VlnPlot(Sobj, features = c("nFeature_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("nCount_RNA"), pt.size = 0.1)
VlnPlot(Sobj, features = c("nCount_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("percent.mt"), pt.size = 0.1)
VlnPlot(Sobj, features = c("percent.mt"), pt.size = 0)

Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message:
“The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
ℹ Please use the `layer` argument instead.
ℹ The deprecated feature was likely used in the Seurat package.
  Please report the issue at <https://github.com/satijalab/seurat/issues>.”
Warning message:
“`PackageCheck()` was deprecated in SeuratObject 5.0.0.
ℹ Please use `rlang::check_installed()` instead.
ℹ The deprecated feature was likely used in the Seurat package.
  Please report the issue at <https://github.com/satijalab/seurat/issues>.”
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image 
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image 
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image 
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image 
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image
No description has been provided for this image
In [10]:
p1 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "percent.mt", raster = F)
p2 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = F)
p1+p2
t_1 <- table(Sobj$orig.ident)
t_1

sample_1 sample_2 sample_3 
   20092    19107    24589
No description has been provided for this image
In [11]:
cutoffs <- data.frame(matrix(NA,nrow = 3, ncol = 3))
colnames(cutoffs) <- c("nF_cutoff", "mt_cutoff", "mt_cutoff_sd")
rownames(cutoffs) <- paste0("sample_", c(1:3))
cell_ind <- list()
for (a in rownames(cutoffs)) {
  cutoffs[a,"nF_cutoff"] <- median(Sobj$nFeature_RNA[Sobj$orig.ident==a]) + sd(Sobj$nFeature_RNA[Sobj$orig.ident==a])*5
  cutoffs[a,"mt_cutoff_sd"] <- median(Sobj$percent.mt[Sobj$orig.ident==a]) + sd(Sobj$percent.mt[Sobj$orig.ident==a])*2
}
In [12]:
cutoffs[,"mt_cutoff"] <- cutoffs[,"mt_cutoff_sd"]
cutoffs[cutoffs[,"mt_cutoff_sd"]>25,"mt_cutoff"] <- 25
cell_ind <- NULL
for (s_name in rownames(cutoffs)) {
  tm_ind <- colnames(Sobj[,Sobj$orig.ident==s_name & Sobj$nFeature_RNA < cutoffs[s_name,"nF_cutoff"] & Sobj$nFeature_RNA>300 & Sobj$percent.mt < cutoffs[s_name,"mt_cutoff"]])
  cell_ind <- c(cell_ind,tm_ind)
}
In [13]:
expressed_cell_count <- data.frame(matrix(NA,nrow = dim(Sobj)[1], ncol = 9))
colnames(expressed_cell_count) <- rownames(cutoffs)
for (i in rownames(cutoffs)) {
  expressed_cell_count[,i] <- rowSums(Sobj@assays$RNA$counts[,Sobj$orig.ident==i]>1)
}
print(table(Sobj$orig.ident))

sample_1 sample_2 sample_3 
   20092    19107    24589
In [14]:
print(colSums(expressed_cell_count>10))
keep_genes <- rowSums(expressed_cell_count>10)>0

sample_1 sample_2 sample_3     <NA>     <NA>     <NA>     <NA>     <NA> 
   12635    12425    12374       NA       NA       NA       NA       NA 
    <NA> 
      NA
In [15]:
Sobj$cell_filtered <- is.element(colnames(Sobj), cell_ind)

for (s_name in rownames(cutoffs)) {
  a0 <- VlnPlot(Sobj[,Sobj$orig.ident==s_name], features = c("nFeature_RNA","nCount_RNA","percent.mt"), pt.size = 0)
  a1 <- FeatureScatter(Sobj[,Sobj$orig.ident==s_name], feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "cell_filtered") + theme(legend.position = "none") + ggtitle(s_name)
  a2<- FeatureScatter(Sobj[,Sobj$orig.ident==s_name], feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "cell_filtered") + theme(legend.position = "none") + ggtitle(s_name)
  a3<- FeatureScatter(Sobj[,Sobj$orig.ident==s_name], feature1 = "nFeature_RNA", feature2 = "percent.mt", group.by = "cell_filtered") + theme(legend.position = "none") + ggtitle(s_name)
  print(a0)
  print(a1+a2+a3)
}

Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image 
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image
No description has been provided for this image 
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [16]:
Sobj <- subset(Sobj, subset = cell_filtered, features = rownames(Sobj)[keep_genes])
dim(Sobj)
  1. 33696
  2. 62860
In [17]:
Sobj[["CITE"]] <- CreateAssayObject(counts = bb[,cell_ind])

Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
In [18]:
t_2 <- table(Sobj$orig.ident)
t_2

sample_1 sample_2 sample_3 
   19869    18861    24130
In [19]:
VlnPlot(Sobj, features = c("nFeature_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("nCount_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("percent.mt"), pt.size = 0)

Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image 
Warning message:
“Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
No description has been provided for this image
No description has been provided for this image
In [20]:
Idents(Sobj) <- Sobj$orig.ident
p1 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "percent.mt", raster = F)
p2 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = F)
p1+p2
No description has been provided for this image
In [21]:
Sobj <- NormalizeData(Sobj, normalization.method = "LogNormalize", scale.factor = 10000)
Sobj <- FindVariableFeatures(Sobj, selection.method = "vst", nfeatures = 2000)
Sobj <- ScaleData(Sobj)
Sobj <- RunPCA(Sobj)
DimPlot(Sobj)
ElbowPlot(Sobj)

Normalizing layer: counts

Finding variable features for layer counts

Centering and scaling data matrix

PC_ 1 
Positive:  Igkc, Iglc3, Ighm, Gm30211, Cr2, Mzb1, Kcnq5, Myc, Mast4, E330020D12Rik 
	   Immp2l, Inpp4b, Maml2, 2010309G21Rik, Iglc1, Ncl, Hs3st1, Dtx1, Arhgap24, Rab30 
	   Plk2, Cacna1e, Txndc5, Mif, Ly6a, Hdac9, Ift140, Rps2, Pik3r4, Iglv3 
Negative:  Fcer1g, Ifitm3, Lyz2, Cd300a, Ifitm2, Csf1r, Alox5ap, Ccl6, Fyb, Clec7a 
	   Ifitm6, S100a4, Metrnl, Clec4a1, Clec4a3, Igsf6, Lgals3, Emilin2, Tnfrsf1a, Gm36161 
	   Sirpa, Gpr141, Cd300c2, Pid1, Gsr, Lrp1, Gm21188, Ifi207, Itgam, Lst1 
PC_ 2 
Positive:  Ace, Pparg, Treml4, Ear2, Cd300e, Cd300ld, Adgre4, S1pr5, Fabp4, Eno3 
	   Pilra, Tgm2, Smpdl3b, Fcgr4, Zfyve9, Slc11a1, Apoc2, Rap1gap2, Pilrb2, Tiam2 
	   Hfe, Arhgef10l, Cd300c2, Trem3, Gcnt2, Havcr2, Pglyrp1, Fgd4, Adrb1, Lilra5 
Negative:  Mki67, Pclaf, Birc5, Kif11, Ccna2, Top2a, Cdca8, Cdk1, Spc24, Tpx2 
	   Cks1b, Cdca3, Racgap1, Stmn1, Kif15, Ube2c, Hmmr, Cenpf, Ccnb2, Prc1 
	   Kif22, H2ac4, Ncapg, H2az1, Clspn, Esco2, Neil3, Bub1, Kif4, Cdca5 
PC_ 3 
Positive:  F13a1, Fn1, Ccr2, AI839979, Ccl9, Cd177, Ly6c2, Stxbp6, Plcb1, Slfn1 
	   C3, Hp, Ly6a2, Wfdc17, Clec4a2, Gas7, Olfm1, Hacd4, Ifi211, Ms4a8a 
	   Msr1, Cd14, S100a6, Thbs1, Mmp8, Mcub, Tgfbi, Mgst1, Arhgef37, S100a4 
Negative:  Pparg, Cd300e, Ace, Treml4, Eno3, Ear2, Tgm2, Cd9, Lair1, Cd300ld 
	   Havcr2, Cd36, S1pr5, Fabp4, Adgre4, Gcnt2, Itgax, Fcgr4, Spn, Hfe 
	   Cdk14, Mki67, Grk3, Pilra, Smpdl3b, Slc12a2, Pclaf, Rap1gap2, Runx2, Ppm1h 
PC_ 4 
Positive:  Mki67, Pclaf, Birc5, Ccna2, Kif11, Cdk1, Spc24, Cdca3, Top2a, Cdca8 
	   Tpx2, Hmmr, Ube2c, Cenpf, Kif15, Ccnb2, Prc1, Kif22, H2ac4, Stmn1 
	   Esco2, Cenpe, Neil3, Clspn, Bub1, Nuf2, Kif4, Rrm2, Ncapg, Cks1b 
Negative:  Cd7, Klrd1, Tcf7, Itk, Skap1, Nkg7, Ms4a4b, Cd3d, Cd3g, Gm2682 
	   Cd3e, Lat, Cd28, Prkcq, Trbc2, Cd247, Il2rb, Ccl5, Cd27, Slc9a9 
	   Il7r, Ccnd1, Klrk1, Srgap3, Prkch, Thy1, Ctsw, Tbc1d4, Tox, Bcl11b 
PC_ 5 
Positive:  Atxn1, Pik3r4, Zc3h12c, Ahr, Myof, Fcrl5, Dtx1, Cr2, Myc, Atf3 
	   Nebl, S1pr3, Plac8, Tank, Ackr3, Odc1, Dph5, Hs3st1, Cd86, Pdia4 
	   Rilpl2, Tbc1d9, Camkmt, AW011738, Dnase1l3, Immp2l, Cybb, Hpgds, Marcks, Jun 
Negative:  Itk, Tcf7, Ms4a4b, Nkg7, Cd3d, Gm2682, Cd3e, Cd3g, Trbc2, Lat 
	   Cd247, Prkcq, Cd28, Cd27, Skap1, Tox, Prkch, Il2rb, Lef1, Thy1 
	   Bcl11b, Trbc1, Themis, Ctsw, Il7r, Camk4, Sh2d1a, Ctla2a, Txk, Ikzf2
No description has been provided for this image
No description has been provided for this image
In [22]:
Sobj <- RunUMAP(Sobj, dims = 1:15)

Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
17:01:14 UMAP embedding parameters a = 0.9922 b = 1.112

17:01:14 Read 62860 rows and found 15 numeric columns

17:01:14 Using Annoy for neighbor search, n_neighbors = 30

17:01:14 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

17:01:22 Writing NN index file to temp file /tmp/RtmpiDhTyy/file5b1f448442d2a

17:01:22 Searching Annoy index using 1 thread, search_k = 3000

17:01:46 Annoy recall = 100%

17:01:47 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

17:01:49 Initializing from normalized Laplacian + noise (using RSpectra)

17:01:55 Commencing optimization for 200 epochs, with 2692262 positive edges

17:01:55 Using rng type: pcg

17:02:36 Optimization finished
In [23]:
DimPlot(Sobj, raster = F)
No description has been provided for this image
In [24]:
library(harmony)

Loading required package: Rcpp
In [25]:
Sobj <- RunHarmony(Sobj,"orig.ident", plot_convergence = T)
harmony_embeddings <- Embeddings(Sobj, "harmony")
harmony_embeddings[1:5, 1:5]

Transposing data matrix

Initializing state using k-means centroids initialization

Warning message:
“Quick-TRANSfer stage steps exceeded maximum (= 3143000)”
Warning message:
“Quick-TRANSfer stage steps exceeded maximum (= 3143000)”
Warning message:
“Quick-TRANSfer stage steps exceeded maximum (= 3143000)”
Warning message:
“Quick-TRANSfer stage steps exceeded maximum (= 3143000)”
Warning message:
“Quick-TRANSfer stage steps exceeded maximum (= 3143000)”
Warning message:
“Quick-TRANSfer stage steps exceeded maximum (= 3143000)”
Warning message:
“Quick-TRANSfer stage steps exceeded maximum (= 3143000)”
Harmony 1/10

Harmony 2/10

Harmony 3/10

Harmony converged after 3 iterations

Warning message:
“The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
ℹ Please use the `layer` argument instead.
ℹ The deprecated feature was likely used in the Seurat package.
  Please report the issue at <https://github.com/satijalab/seurat/issues>.”
A matrix: 5 × 5 of type dbl
harmony_1harmony_2harmony_3harmony_4harmony_5
sample_1_AAACCAAAGACAATTC-11.449980-1.9607286-2.1336740-0.1322129 3.99104252
sample_1_AAACCAAAGACGACGA-12.540081 0.9835546 1.1427124 0.3662081-0.06780086
sample_1_AAACCAAAGATGTATC-12.281648 1.0801510-0.4316720-1.0081381 2.18528764
sample_1_AAACCAAAGCCAACCG-13.344910 1.1810851 0.1232371 1.0964129-0.49868922
sample_1_AAACCAAAGCCGATTA-11.612576 0.1193344 0.3849787-0.7293387 1.97029635
No description has been provided for this image
In [26]:
DimPlot(Sobj, reduction = "harmony", group.by = "orig.ident", raster = F)
No description has been provided for this image
In [27]:
Sobj <- RunUMAP(Sobj, reduction = "harmony", dims = 1:15)
Sobj <- FindNeighbors(Sobj, reduction = "harmony", dims = 1:15)
Sobj <- FindClusters(Sobj, resolution = 0.5)

17:03:56 UMAP embedding parameters a = 0.9922 b = 1.112

17:03:56 Read 62860 rows and found 15 numeric columns

17:03:56 Using Annoy for neighbor search, n_neighbors = 30

17:03:56 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

17:04:03 Writing NN index file to temp file /tmp/RtmpiDhTyy/file5b1f4286fb1ac

17:04:03 Searching Annoy index using 1 thread, search_k = 3000

17:04:28 Annoy recall = 100%

17:04:29 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

17:04:31 Initializing from normalized Laplacian + noise (using RSpectra)

17:04:36 Commencing optimization for 200 epochs, with 2705776 positive edges

17:04:36 Using rng type: pcg

17:05:17 Optimization finished

Computing nearest neighbor graph

Computing SNN


Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 62860
Number of edges: 1625800

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8656
Number of communities: 17
Elapsed time: 21 seconds
In [28]:
DimPlot(Sobj, raster = F)
DimPlot(Sobj, raster = F, group.by = "orig.ident")
No description has been provided for this image
No description has been provided for this image
In [29]:
Sobj <- ScaleData(Sobj, features = rownames(Sobj))

Centering and scaling data matrix

Warning message:
“Different features in new layer data than already exists for scale.data”
In [30]:
source("~/projects/MK_Rcode/setZ.R")
ABC_marker_mouse <- c("Cd19", "Cd80", "Cd86", "Fas", "Fcrla", "Fcrlb", "Fgl2", "Cxcl10", "Cxcr3", "Nkg7", "Itgam", "Itgb2l", "Itgb2", "Tbx21", "Zbtb32", "Fcer2a", "Cr2")
length(ABC_marker_mouse)
is.element(ABC_marker_mouse,rownames(Sobj))
ABC_marker_mouse <- ABC_marker_mouse[is.element(ABC_marker_mouse,rownames(Sobj))]
length(ABC_marker_mouse)
17
  1. TRUE
  2. TRUE
  3. TRUE
  4. TRUE
  5. TRUE
  6. TRUE
  7. TRUE
  8. TRUE
  9. TRUE
  10. TRUE
  11. TRUE
  12. TRUE
  13. TRUE
  14. TRUE
  15. TRUE
  16. TRUE
  17. TRUE
17
In [31]:
Sobj[["ABC_score"]] <- setZ(ABC_marker_mouse,rownames(Sobj@assays$RNA$scale.data),Sobj@assays$RNA$scale.data)
In [32]:
LSP1_gene_mouse <- c("Cebpa", "Cebpb", "Il1b", "Ccl3", "Tnf", "Ccl4", "Ccl5", "Spp1", "Cd7", "Msr1", "Itgam", "Kit", "Cd14", "Csf3r", "Il1r2", "Csf1r", "Itga1", "Itga5", "Itgb5", "Itgb1", "Camp", "Thbs1", "C3", "Mefv", "Fos", "Cybb", "Card9", "Irf7", "Oas3", "Nlrp3", "Ctsb", "Ikbke", "Mpo")
length(LSP1_gene_mouse)
33
In [33]:
LSP1_gene_sub <- LSP1_gene_mouse[is.element(LSP1_gene_mouse,rownames(Sobj))]
length(LSP1_gene_sub)
33
In [34]:
Sobj[["LSP1_score"]] <- setZ(LSP1_gene_sub,rownames(Sobj@assays$RNA$scale.data),Sobj@assays$RNA$scale.data)
In [35]:
Sobj$LSP1_score_pos <- Sobj$LSP1_score > 0
Sobj$ABC_score_pos <- Sobj$ABC_score > 0
FeaturePlot(Sobj, features = "LSP1_score") + ggtitle("LSP1 score")
FeaturePlot(Sobj, features = "ABC_score") + ggtitle("ABC score")
p0 <- DimPlot(Sobj, group.by = "LSP1_score_pos") + ggtitle("LRB")
p0
p0_1 <- DimPlot(Sobj, group.by = "ABC_score_pos") + ggtitle("ABC")
p0_1
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [36]:
if (!file.exists(paste0(path_out,"Lee_WT_harmony_052025.rds"))) {
    saveRDS(Sobj, paste0(path_out,"Lee_WT_harmony_052025.rds"))
    saveRDS(Sobj@meta.data, paste0(path_out,"Lee_WT_metadata_052025.rds"))
}
In [37]:
marker_list <- c("Myc", "Ebf1", "Pax5", "Dtx1", "Bank1", "Cd55", "Cd19", "Cd93", "Cd24a", "Cd38", "Ighm", "Ighd", "Ighg", "Cd23", "Cd21", "Sell", "Cd5", "Cd43", "Cd27", "Cd80", "Cd73", "Pdcd1lg2", "Bcl6", "Aicda", "Mki67", "Cxcr4", "Cd83", "Cd86", "Cd74", "Prdm1", "Xbp1", "Irf4", "Cd138", "Jchain", "Mzb1")
In [38]:
Idents(Sobj) <- Sobj$seurat_clusters
pdf(paste0(path_out,"marker_dotplot_052125.pdf"), width = 20, height = 6)
DotPlot(Sobj, features = marker_list) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)
dev.off()

Warning message:
“The following requested variables were not found: Ighg, Cd23, Cd21, Cd43, Cd73, Cd138”
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
pdf: 2
In [39]:
if (!file.exists(paste0(path_out,"marker_list_WT_052125.rds"))) {
    all.markers <- FindAllMarkers(Sobj, test.use = "MAST", min.pct = 0.1, logfc.threshold = 0.58, verbose = F, only.pos = T)
    saveRDS(all.markers, file = paste0(path_out,"marker_list_WT_052125.rds"))
}
if (!exists("all.markers")) {
    all.markers <- readRDS(paste0(path_out,"marker_list_WT_052125.rds"))
}
In [40]:
marker_200 <- all.markers %>%
    group_by(cluster) %>%
    filter(avg_log2FC > 0.58) %>%
    slice_head(n = 200) %>%
    ungroup()
In [41]:
write.csv(marker_200, paste0(path_out, "top200_marker_genes.csv"))
In [42]:
Idents(Sobj) <- Sobj$seurat_clusters
Sobj <- RenameIdents(Sobj,
                     "0" = "Naive Follicular B", "1" = "Naive Follicular B", "2" = "Naive Follicular B", "3" = "Naive Follicular B",
                     "4" = "Transitional B (T2)", "5" = "Naive Follicular B", "6" = "Late transitional B", "7" = "Naive Follicular B",
                     "8" = "Plasma cell", "9" = "Myeloid-like GC B", "10" = "Myeloid-like GC B", "11" = "Myeloid-like GC B",
                     "12" = "Mki67+ GC DZ B", "13" = "T/B Doublet", "14" = "GC LZ B", "15" = "Naive Follicular B", "16" = "Plasmablast"
                     )
Sobj$annotation <- Idents(Sobj)
In [43]:
d1 <- DimPlot(Sobj, group.by = "seurat_clusters", label = T, repel = T)
d2 <- DimPlot(Sobj, group.by = "annotation", label = T, repel = T)
d1
d2

pdf(paste0(path_out,"umap_052125.pdf"), width = 6, height = 4)
 d1 + ggtitle("")
 d2 + ggtitle("")
dev.off()
pdf(paste0(path_out,"umap_053025v1.pdf"), width = 4, height = 4)
d1 + ggtitle("") + NoLegend()
d2 + ggtitle("") + NoLegend()
DimPlot(Sobj, group.by = "annotation", label = F) + ggtitle("") + NoLegend()
p0
p0 + ggtitle("") + NoLegend()
p0_1
p0_1 + ggtitle("") + NoLegend()
dev.off()
No description has been provided for this image
pdf: 2
pdf: 2
No description has been provided for this image
In [44]:
RidgePlot(Sobj, features = c("LSP1_score", "ABC_score"), ncol = 2)
d3 <- RidgePlot(Sobj, features = c("LSP1_score")) + NoLegend() + ggtitle("")
d4 <- RidgePlot(Sobj, features = c("ABC_score")) + NoLegend() + ggtitle("")

Picking joint bandwidth of 0.385

Picking joint bandwidth of 0.355
No description has been provided for this image
In [45]:
pdf(paste0(path_out, "RidgePlot_LSP1_ABC_052225.pdf"), width = 4, height = 4)
d3
d4
dev.off()

Picking joint bandwidth of 0.385

Picking joint bandwidth of 0.355
pdf: 2
In [46]:
celltype.prop <- as.data.frame(table(Sobj$annotation))
celltype.prop <- celltype.prop %>%
    rename(celltype = Var1) %>%
    mutate(Percent = Freq / sum(Freq)*100)
d5 <- ggplot(celltype.prop, aes(x = "", y = Percent, fill = celltype))+
    geom_bar(stat = "identity", width = 0.9) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), axis.title.x = element_blank(), plot.margin = margin(10, 10, 20, 10)) + 
    NoLegend()
d5
No description has been provided for this image
In [47]:
pdf(paste0(path_out,"stacked_bar_WT_052225.pdf"), width = 1.5, height = 4)
d5
dev.off()
pdf: 2
In [48]:
celltype.prop
A data.frame: 9 × 3
celltypeFreqPercent
<fct><int><dbl>
Naive Follicular B 4906378.0512249
Transitional B (T2) 5054 8.0400891
Late transitional B 2659 4.2300350
Plasma cell 1801 2.8650970
Myeloid-like GC B 2930 4.6611518
Mki67+ GC DZ B 488 0.7763283
T/B Doublet 413 0.6570156
GC LZ B 342 0.5440662
Plasmablast 110 0.1749920
In [49]:
table(Sobj$LSP1_score_pos)
table(Sobj$ABC_score_pos)

FALSE  TRUE 
47330 15530 
FALSE  TRUE 
34780 28080
In [50]:
LSP1_pos_ind <- Sobj$LSP1_score_pos==TRUE & Sobj$ABC_score_pos==FALSE
ABC_pos_ind <- Sobj$LSP1_score_pos==FALSE & Sobj$ABC_score_pos==TRUE
DP_ind <- Sobj$LSP1_score_pos==TRUE & Sobj$ABC_score_pos==TRUE
DN_ind <- Sobj$LSP1_score_pos==FALSE & Sobj$ABC_score_pos==FALSE

Sobj$anno_ABC_LSP1 <- NA
Sobj$anno_ABC_LSP1[LSP1_pos_ind] <- "LRB"
Sobj$anno_ABC_LSP1[ABC_pos_ind] <- "ABC"
Sobj$anno_ABC_LSP1[DP_ind] <- "LRB & ABC"
Sobj$anno_ABC_LSP1[DN_ind] <- "Classical B"
unique(Sobj$anno_ABC_LSP1)
Idents(Sobj) <- factor(Sobj$anno_ABC_LSP1, levels = c("LRB", "ABC", "LRB & ABC", "Classical B"))
Sobj$anno_ABC_LSP1 <- Idents(Sobj)
levels(Idents(Sobj))
  1. 'Classical B'
  2. 'ABC'
  3. 'LRB'
  4. 'LRB & ABC'
  1. 'LRB'
  2. 'ABC'
  3. 'LRB & ABC'
  4. 'Classical B'
In [51]:
options(repr.plot.width=10, repr.plot.height=8)
Idents(Sobj) <- Sobj$anno_ABC_LSP1
levels(Sobj$anno_ABC_LSP1)
Idents(Sobj) <- Sobj$anno_ABC_LSP1
my_colors <- c("LRB" = "#ff0000", "ABC" = "#ffff00", "LRB & ABC" = "#00b0ff", "Classical B" = "#BFBFBF")
d6 <- DimPlot(Sobj, reduction = "umap", cols = my_colors)
d6
  1. 'LRB'
  2. 'ABC'
  3. 'LRB & ABC'
  4. 'Classical B'
No description has been provided for this image
In [52]:
if (!file.exists(paste0(path_out,"Lee_WT_harmony_080425.rds"))) {
    saveRDS(Sobj, paste0(path_out,"Lee_WT_harmony_080425.rds"))
}
In [53]:
library(ggpointdensity)
library(scales)
In [54]:
df <- FetchData(Sobj, vars = c("umap_1", "umap_2", "LSP1_score_pos", "ABC_score_pos", "anno_ABC_LSP1"))
df_LSP1 <- df %>% filter(LSP1_score_pos)
df_ABC <- df %>% filter(ABC_score_pos)
In [55]:
options(repr.plot.width=6, repr.plot.height=5)
t1 <- ggplot(df_LSP1, aes(x = umap_1, y = umap_2)) +
    geom_pointdensity(method = "neighbors", size = 0.1) +
      scale_color_gradientn(
        colors = c("grey", "grey", "red"),
        values = rescale(c(0, 350, 1200)),
        limits = c(0, 1200)
      ) +
    theme_classic()

t2 <- ggplot(df_ABC, aes(x = umap_1, y = umap_2)) +
  geom_pointdensity(method = "neighbors", size = 0.1) +
  scale_color_gradientn(
    colors = c("grey", "grey", "red"),
    values = rescale(c(0, 350, 1200)),
    limits = c(0, 1200)
  ) +
    theme_classic()
t1
t2
No description has been provided for this image
No description has been provided for this image
In [56]:
pdf(paste0(path_out,"LRB_ABC_umap_density_080425.pdf"), width = 6, height = 5)
t1
t2
dev.off()
pdf: 2
In [57]:
t1_1 <- t1 + 
labs(x = "", y = "") + 
theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(),
    axis.ticks.y = element_line()
  )
t2_1 <- t2 + 
labs(x = "", y = "") + 
theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(),
    axis.ticks.y = element_line()
  )
t1_1
t2_1
No description has been provided for this image
No description has been provided for this image
In [70]:
tiff(filename = paste0(path_out,"LRB_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(t1_1)
dev.off()

tiff(filename = paste0(path_out,"ABC_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(t2_1)
dev.off()
pdf: 2
pdf: 2
In [69]:
path_out
'~/projects/2025/LeeNaeun/out/'
In [59]:
library(clusterProfiler)
library(enrichplot)
library(msigdbr)
library(stringr)
set.seed(2025)
hallmark_geneset <- msigdbr(species = "mouse",  category = "H")


clusterProfiler v4.14.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

G Yu. Thirteen years of clusterProfiler. The Innovation. 2024,
5(6):100722


Attaching package: ‘clusterProfiler’


The following object is masked from ‘package:stats’:

    filter


enrichplot v1.26.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

Guangchuang Yu. Using meshes for MeSH term enrichment and semantic
analyses. Bioinformatics. 2018, 34(21):3766-3767,
doi:10.1093/bioinformatics/bty410

Using human MSigDB with ortholog mapping to mouse. Use `db_species = "MM"` for mouse-native gene sets.
This message is displayed once per session.
Warning message:
“The `category` argument of `msigdbr()` is deprecated as of msigdbr 10.0.0.
ℹ Please use the `collection` argument instead.”
In [60]:
hallmark_geneset$gs_name_short <- str_extract(hallmark_geneset$gs_name, "(?<=HALLMARK_).*$")
head(hallmark_geneset)
A tibble: 6 × 27
gene_symbolncbi_geneensembl_genedb_gene_symboldb_ncbi_genedb_ensembl_genesource_genegs_idgs_namegs_collection⋯gs_urldb_versiondb_target_speciesortholog_taxon_idortholog_sourcesnum_ortholog_sourcesentrez_genegs_catgs_subcatgs_name_short
<chr><chr><chr><chr><chr><chr><chr><chr><chr><chr>⋯<chr><chr><chr><int><chr><dbl><chr><chr><chr><chr>
Abca111303ENSMUSG00000015243ABCA119 ENSG00000165029ABCA1M5905HALLMARK_ADIPOGENESISH⋯2025.1.HsHS10090EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|Treefam 1111303HADIPOGENESIS
Abcb874610ENSMUSG00000028973ABCB811194ENSG00000197150ABCB8M5905HALLMARK_ADIPOGENESISH⋯2025.1.HsHS10090EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|PhylomeDB|Treefam1274610HADIPOGENESIS
Acaa252538ENSMUSG00000036880ACAA210449ENSG00000167315ACAA2M5905HALLMARK_ADIPOGENESISH⋯2025.1.HsHS10090EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|Treefam 1152538HADIPOGENESIS
Acadl11363ENSMUSG00000026003ACADL33 ENSG00000115361ACADLM5905HALLMARK_ADIPOGENESISH⋯2025.1.HsHS10090EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|PhylomeDB|Treefam1211363HADIPOGENESIS
Acadm11364ENSMUSG00000062908ACADM34 ENSG00000117054ACADMM5905HALLMARK_ADIPOGENESISH⋯2025.1.HsHS10090EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|PhylomeDB|Treefam1211364HADIPOGENESIS
Acads11409ENSMUSG00000029545ACADS35 ENSG00000122971ACADSM5905HALLMARK_ADIPOGENESISH⋯2025.1.HsHS10090EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|Treefam 1111409HADIPOGENESIS
In [63]:
if (!file.exists(paste0(path_out,"LSP1_vs_ABC.marker.rds"))) {
    LSP1_vs_ABC.marker <- FindMarkers(Sobj, ident.1 = "LRB", ident.2 = "ABC", min.pct = 0.1, logfc.threshold = 0.1, test.use = "MAST")
    saveRDS(LSP1_vs_ABC.marker, file = paste0(path_out,"LSP1_vs_ABC.marker.rds"))
}

if (!exists("LSP1_vs_ABC.marker")) {
    LSP1_vs_ABC.marker <- readRDS(paste0(path_out,"LSP1_vs_ABC.marker.rds"))
}
In [64]:
sort_glist <- LSP1_vs_ABC.marker$avg_log2FC
names(sort_glist) <- rownames(LSP1_vs_ABC.marker)
sort_glist <- sort(sort_glist, decreasing = T)
gsea_result.LSP1vsABS <- GSEA(geneList = sort_glist, TERM2GENE = select(hallmark_geneset, gs_name_short, gene_symbol), seed = T, pvalueCutoff = 0.99)

using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).


preparing geneSet collections...

GSEA analysis...

leading edge analysis...

done...
In [65]:
aa <- dotplot(gsea_result.LSP1vsABS, showCategory=10, split=".sign") + facet_grid(.~.sign)
aa
saveRDS(gsea_result.LSP1vsABS, file = paste0(path,"gsea_.rds"))
No description has been provided for this image
In [66]:
head(gsea_result.LSP1vsABS@result)
write.csv(gsea_result.LSP1vsABS@result, file = paste0(path_out,"gsea_mouse_LRB_vs_ABS.csv"))
A data.frame: 6 × 11
IDDescriptionsetSizeenrichmentScoreNESpvaluep.adjustqvaluerankleading_edgecore_enrichment
<chr><chr><int><dbl><dbl><dbl><dbl><dbl><dbl><chr><chr>
INTERFERON_GAMMA_RESPONSEINTERFERON_GAMMA_RESPONSEINTERFERON_GAMMA_RESPONSE750.67685372.0741812.107233e-050.00074924220.0004907317551tags=51%, list=18%, signal=43%Ifitm3/Ifitm2/Irf7/Bst2/Isg15/Pim1/Zbp1/Oasl1/Xaf1/Rnf213/Isg20/Eif2ak2/Nod1/Nfkbia/Herc6/Tdrd7/Trafd1/Mvp/Eif4e3/Parp14/Stat1/Socs3/Icam1/Nlrc5/Irf9/Epsti1/Myd88/Tnfaip3/Ifi35/Il10ra/Ifnar2/Psme2/Slamf7/Casp4/Ogfr/Znfx1/Casp1/Irf1
TNFA_SIGNALING_VIA_NFKBTNFA_SIGNALING_VIA_NFKB TNFA_SIGNALING_VIA_NFKB 910.64405681.9922523.329965e-050.00074924220.0004907317746tags=66%, list=24%, signal=51%Tnf/Ccl4/Ptpre/Cebpb/Atf3/Klf10/Cd44/Ldlr/Sgk1/Sat1/Trib1/Dusp1/Egr1/Marcks/Nfkbia/Nr4a3/Egr2/Ier2/Fos/Traf1/Zfp36/Socs3/Pdlim5/Icam1/Klf4/Fosb/Rnf19b/Btg2/Dusp5/Tnfaip3/Jun/Nr4a2/Nr4a1/Maff/Mxd1/Pfkfb3/Nfe2l2/Kdm6b/Map2k3/Gadd45b/Map3k8/Atp2b1/Junb/Hes1/Irf1/Plek/Ptger4/B4galt5/Ppp1r15a/Ninj1/Sqstm1/Il6st/Mcl1/Tsc22d1/Smad3/Rhob/Stat5a/Dennd5a/Myc/Tank
INTERFERON_ALPHA_RESPONSEINTERFERON_ALPHA_RESPONSEINTERFERON_ALPHA_RESPONSE400.73555142.1087759.441205e-050.00141618070.0009275570342tags=40%, list=11%, signal=36%Ifitm3/Ifitm2/Irf7/Bst2/Isg15/Oasl1/Isg20/Eif2ak2/Tent5a/Herc6/Tdrd7/Trafd1/Parp14/Irf9/Parp9/Epsti1
IL6_JAK_STAT3_SIGNALINGIL6_JAK_STAT3_SIGNALING IL6_JAK_STAT3_SIGNALING 360.72927692.0685791.512185e-040.00143676350.0009410381375tags=47%, list=12%, signal=42%Tnf/Il1r2/Csf2ra/Cd9/Pim1/Cd44/Cd36/Tnfrsf1b/Ifngr1/Ebi3/Il17ra/Pik3r5/Stat1/Socs3/Irf9/Myd88/Jun
APOPTOSISAPOPTOSIS APOPTOSIS 500.68136052.0055671.915685e-040.00143676350.0009410381551tags=48%, list=18%, signal=40%Ifitm3/Tnf/Lgals3/Atf3/Cd44/H1f0/Gpx1/Sat1/Isg20/Ifngr1/Pmaip1/Pea15a/Ank/Btg2/Jun/Dap/Sod1/Gadd45b/Tspo/Casp4/Ppt1/Rara/Casp1/Irf1
COMPLEMENTCOMPLEMENT COMPLEMENT 580.66706811.9835641.667369e-040.00143676350.0009410381120tags=21%, list=4%, signal=20% Fcer1g/Lgals3/Irf7/Ctsb/Msrb1/Cebpb/Anxa5/Pim1/Dgkh/Gngt2/Atox1/Cd36
In [67]:
path_out
'~/projects/2025/LeeNaeun/out/'
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [68]:
sessionInfo()

R version 4.4.2 (2024-10-31)
Platform: x86_64-conda-linux-gnu
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS/LAPACK: /home/vincent/aprojects/Yeonjoo/anaconda3/envs/r_env_restored/lib/libopenblasp-r0.3.29.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=ko_KR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=ko_KR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=ko_KR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=ko_KR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Los_Angeles
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] stringr_1.5.1          msigdbr_25.1.0         enrichplot_1.26.6     
 [4] clusterProfiler_4.14.6 scales_1.4.0           ggpointdensity_0.2.0  
 [7] future_1.58.0          harmony_1.2.3          Rcpp_1.1.0            
[10] dplyr_1.1.4            ggplot2_3.5.2          Seurat_5.3.0          
[13] SeuratObject_5.1.0     sp_2.2-0              

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22        splines_4.4.2           later_1.4.2            
  [4] pbdZMQ_0.3-14           ggplotify_0.1.2         tibble_3.3.0           
  [7] R.oo_1.27.1             polyclip_1.10-7         fastDummies_1.7.5      
 [10] lifecycle_1.0.4         globals_0.18.0          lattice_0.22-7         
 [13] MASS_7.3-65             magrittr_2.0.3          plotly_4.11.0          
 [16] ggtangle_0.0.7          httpuv_1.6.16           sctransform_0.4.2      
 [19] spam_2.11-1             spatstat.sparse_3.1-0   reticulate_1.42.0      
 [22] cowplot_1.1.3           pbapply_1.7-2           DBI_1.2.3              
 [25] RColorBrewer_1.1-3      abind_1.4-8             zlibbioc_1.52.0        
 [28] Rtsne_0.17              purrr_1.0.4             R.utils_2.13.0         
 [31] BiocGenerics_0.52.0     yulab.utils_0.2.0       GenomeInfoDbData_1.2.13
 [34] IRanges_2.40.1          S4Vectors_0.44.0        ggrepel_0.9.6          
 [37] irlba_2.3.5.1           listenv_0.9.1           spatstat.utils_3.1-4   
 [40] tidytree_0.4.6          goftest_1.2-3           RSpectra_0.16-2        
 [43] spatstat.random_3.4-1   fitdistrplus_1.2-4      parallelly_1.45.0      
 [46] codetools_0.2-20        DOSE_4.0.1              tidyselect_1.2.1       
 [49] aplot_0.2.8             UCSC.utils_1.2.0        farver_2.1.2           
 [52] matrixStats_1.5.0       stats4_4.4.2            base64enc_0.1-3        
 [55] spatstat.explore_3.4-3  jsonlite_2.0.0          progressr_0.15.1       
 [58] ggridges_0.5.6          survival_3.8-3          systemfonts_1.2.3      
 [61] tools_4.4.2             treeio_1.30.0           ragg_1.3.3             
 [64] ica_1.0-3               glue_1.8.0              gridExtra_2.3          
 [67] qvalue_2.38.0           GenomeInfoDb_1.42.3     IRdisplay_1.1          
 [70] withr_3.0.2             fastmap_1.2.0           digest_0.6.37          
 [73] gridGraphics_0.5-1      R6_2.6.1                mime_0.13              
 [76] textshaping_1.0.1       scattermore_1.2         GO.db_3.20.0           
 [79] Cairo_1.6-2             tensor_1.5.1            spatstat.data_3.1-6    
 [82] RSQLite_2.4.1           R.methodsS3_1.8.2       RhpcBLASctl_0.23-42    
 [85] tidyr_1.3.1             generics_0.1.4          data.table_1.17.6      
 [88] httr_1.4.7              htmlwidgets_1.6.4       uwot_0.2.3             
 [91] pkgconfig_2.0.3         gtable_0.3.6            blob_1.2.4             
 [94] lmtest_0.9-40           XVector_0.46.0          htmltools_0.5.8.1      
 [97] fgsea_1.32.4            dotCall64_1.2           Biobase_2.66.0         
[100] png_0.1-8               spatstat.univar_3.1-3   ggfun_0.1.9            
[103] reshape2_1.4.4          uuid_1.2-1              curl_6.0.1             
[106] nlme_3.1-168            repr_1.1.7              zoo_1.8-14             
[109] cachem_1.1.0            KernSmooth_2.23-26      parallel_4.4.2         
[112] miniUI_0.1.2            vipor_0.4.7             AnnotationDbi_1.68.0   
[115] ggrastr_1.0.2           pillar_1.11.0           grid_4.4.2             
[118] vctrs_0.6.5             RANN_2.6.2              promises_1.3.3         
[121] xtable_1.8-4            cluster_2.1.8.1         beeswarm_0.4.0         
[124] evaluate_1.0.4          cli_3.6.5               compiler_4.4.2         
[127] rlang_1.1.6             crayon_1.5.3            future.apply_1.20.0    
[130] labeling_0.4.3          plyr_1.8.9              fs_1.6.6               
[133] ggbeeswarm_0.7.2        stringi_1.8.7           viridisLite_0.4.2      
[136] deldir_2.0-4            BiocParallel_1.40.2     babelgene_22.9         
[139] assertthat_0.2.1        Biostrings_2.74.1       lazyeval_0.2.2         
[142] spatstat.geom_3.4-1     GOSemSim_2.32.0         Matrix_1.7-3           
[145] IRkernel_1.3.2          RcppHNSW_0.6.0          patchwork_1.3.1        
[148] bit64_4.6.0-1           KEGGREST_1.46.0         shiny_1.11.1           
[151] ROCR_1.0-11             igraph_2.0.3            memoise_2.0.1          
[154] ggtree_3.14.0           fastmatch_1.1-6         bit_4.6.0              
[157] gson_0.1.0              ape_5.8-1
In [ ]: